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The multiconfigurational time-dependent Hartree-Fock method (MCTDHF) is applied for simula- 
tions of the two-photon ionization of Helium. We present results for the single- and double ionization 
from the groundstate for photon energies in the non-sequential regime, and compare them to di- 
rect solutions of the Schrodinger equation using the time-dependent (full) Configuration Interaction 
method (TDCI). We find that the single-ionization is accurately reproduced by MCTDHF, whereas 
the double ionization results correctly capture the main trends of TDCI. 



I. INTRODUCTION 

The photoionization of Helium is among the best stud- 
ied processes involving correlated systems in external 
laser fields. Its investigation dates back over 40 years PQ 
and particularly the last two decades have seen a lot of 
clarifying work on the experimental as well as on the the- 
oretical front. Owing to these efforts the single-photon 
ionization is nowadays considered as well understood [5] . 
Several experiments have been performed, e.g. J3[ H], 
theoretical works provided cross sections which show an 
excellent agreement with the measured data, e.g. [5], 
and also time-resolved theoretical results are becoming 
available [B]. 

In recent time, the two-photon ionization has gained a 
lot of attention, being the simplest multiphoton process 
in which correlation effects play a significant role. De- 
spite the large number of theoretical investigations, see 
e.g. [5[ l7llll|. it is still far from being understood, which 
is reflected by several contradicting predictions for differ- 
ential and total cross sections. Only a few experiments 
have been performed up to now [T^l - TH] . the results of 
which were not yet able to definitely decide on the ori- 
gins of the theoretical discrepancies. 

In this work we present new theoretical results on the 
two-photon ionization at photon energies in the direct 
regime, which ranges roughly in between 39.5 eV - half 
of the total groundstate energy of Helium - and 54.4 eV 
- the second ionization potential. In this regime, double 
ionization can occur only if the electrons share the en- 
ergy of the two photons via their Coulomb interaction. 
The sequential double ionization process, in which the 
electrons are ionized successively, separated at least by 
the relaxation-time of the atom, is not possible with only 
two photons. 

The simulation of Helium in external laser fields is 
most directly performed through the time-dependent 
Schrodinger equation (TDSE). The theoretical tools for 
a direct solution of the TDSE are manifold; most of them 
employ an expansion of the two-particle wavefunction in 
terms of coupled spherical harmonics and arrive at a set 
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of equations for the radial functions. The latter is solved 
either on a two-dimensional grid [12] or by expansion in a 
basis set like Sturmian functions [5], B-splines [8], or the 
discrete variable representation [TT]. On the other hand, 
time-independent approaches have been successfully at- 
tempted, for instance the R-matrix Floquet method |16) . 
and lowest-order perturbation theory (LOPT) [10l fTTj . 
Last but not least, approximate methods have been ap- 
plied, such as the time-dependent density functional the- 
ory (TDDFT) pS [19] or time-dependent Hartree-Fock 
(TDHF) [2Ql El]- They have the advantage of being ap- 
plicable to larger atoms, but suffer the main drawback 
that their accuracy is difficult to estimate. 

In this work, we apply the multiconfigurational time- 
dependent Hartree-Fock (MCTDHF) method, which in 
principle provides a control of the accuracy. As the name 
suggests, MCTDHF can be considered as an extension of 
the time-dependent Hartree-Fock scheme, in which the 
wavefunction is expressed by several configurations in- 
stead of a single one. The origin of the method dates 
back to the 1980's, where time-dependent multiconfigu- 
rational schemes were considered for the first time [221 - 
125] . A particularly successful approach is the multicon- 
figurational time-dependent Hartree (MCTDH) method 
introduced by Meyer, Manthe, Cederbaum and coworkers 
[2"SlI28| . This method has later been naturally extended 
to MCTDHF, where the antisymmetry of the wavefunc- 
tion is explicitly enforced by using Slater determinants 
instead of Hartree products as basis states [29H31] . For a 
recent overview on MCTDH and MCTDHF, see the book 

m 

The formulation used in the present work resembles 
the one given by Alon et al. |33j . who derived equations 
based on the second quantization formalism, which, as 
an aside, have also successfully been applied to the treat- 
ment of bosonic particles |34j . Here, we present a spin- 
restricted version of the fermionic MCTDHF equations, 
which is based upon the same principles, but formulated 
in terms of spatial orbitals rather than spin-orbitals. 

In connection with photoionization, MCTDHF has re- 
cently been applied to molecular hydrogen [351 136| resp. 
to a version thereof with a screened interaction po- 
tential [37]. Some works have also investigated mod- 
els of reduced dimensionality, see e.g. Refs. [3"8rl4"U] 



2 



for one-dimensional and Ref. [H] for two-dimensional 
model systems. However, to the best of our knowledge, 
MCTDHF has not been applied to a treatment of the 
three-dimensional Helium atom in external laser fields 
so far. In this work we perform such an investigation, 
the focus of which lies mainly on a critical assessment 
of the capabilities of the method for the description of 
correlated photoionization processes. Therefore, we com- 
pare the MCTDHF results with the time-dependent (full) 
Configuration Interaction results, which mark the fully 
correlated and thus best possible results for a specified 
single-particle basis. 

The paper is organized as follows: In section II, we give 
a detailed overview on the MCTDHF method in general, 
our chosen discretization schemes and the numerical im- 
plementation. In section III, we present the numerical 
results for the groundstate and for the two-photon ion- 
ization of Helium. Thereafter, in section IV, we give 
a discussion and outlook. Finally, in the appendix we 
present the derivation of the spin-restricted MCTDHF 
equations. 



II. THEORY 

A. General considerations 

In this work we are concerned with the solution of the 
time-dependent electronic Schrodinger equation 



ifltl*) =H(t)\9), 
with a Hamiltonian (in atomic units) given by 

N 
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It describes the motion of N electrons in the field 
of an atomic core with charge Z in the static Born- 
Oppenheimer approximation. The electrons interact 
with each other via the two-body Coulomb repulsion, 
and the operator V ex t(t) incorporates the action of an 
external electromagnetic field. In second quantization 
the Hamiltonian reads H21 
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with the single- and double excitation operators 
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set up as linear combinations of the spin-orbital excita- 
tion operators dp}, which destroy (create) a particle in 



the time- dependent state | 4> p<y {t) ). The Hamiltonian (13J) 
implies a spin-restricted treatment, in which the orbitals 
are independent of their respective spin-projection, i.e. 
I 4>pa{t) ) = | 4>pp{t) ) —'■ I 4>p(t) )• The matrix elements of 
the one-body Hamiltonian and the two-body Coulomb 
operator in this basis are time-dependent as well and 
given by 

h pq (t) = [dr ^(r) A - - + V ext (t)\ 0,(r) , (6) 
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where we indicate only the explicit time-dependence due 
to the external field. The interaction V ex t(t) with the 
external field £{t) is considered in dipole approximation 
in the length gauge 



Kxt(i) = £{t) 



(8) 



Another common choice is the velocity gauge. The per- 
formance of the two gauges is well discussed in the liter- 
ature, see e.g. [131131], and will not be analyzed here. 



B. The MCTDHF equations 

In the MCTDHF method, the wavefunction is ex- 
pressed by a linear superposition of time- dependent 
Slater determinants, 

| *) = ^Cnft) I ni a) n\p, • • • , n Ma , n M p; t) , (9) 
ii 

which are written in occupation number representation 
with respect to the spin-orbitals {| 4> pa {t) )}i< p <m, from 
which they directly inherit their time-dependence. The 
sum runs over the complete set of (jf)( N ) Slater de- 
terminants, that can be constructed from N a particles 
with spin-projection a, Np particles with spin-projection 
j3 and M spatial orbitals. Their spin-projection is conse- 
quently restricted to S z — (N a — Np)/2. Alternatively, 
in Eq. (9) we could also use an expansion in configura- 
tion state functions, which are not only eigenfunctions 
of the spin-projection operator (as Slater determinants 
are) , but also of the total spin operator [12] . 

The MCTDHF ansatz is characterized by two sets of 
time-dependent variational parameters, namely by the 
expansion coefficients C n (t) and by the orbitals | 4> p (t) ), 
which determine the actual many-body basis. Their 
equations of motion are found with the Lagrange formu- 
lation of the time-dependent variational principle which 
requires a minimization of the action functional [331 134j 
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with respect to the variational parameters. The Lagrange 
multipliers (ikl(t) are introduced to ensure the orthonor- 
mality of the orbitals during the temporal evolution. The 
procedure of taking functional derivatives is performed 
in detail in Appendix [X] It results in the spin-restricted 
MCTDHF equations 



iC a (t) = ^(n|i?(t)|m)C m (t) : 



^ pqrs 
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a set of coupled first-order differential equation for the 
temporal evolution of the variational parameters. Equa- 
tion (111, here denoted as wavefunction equation, is 
just the Schrodinger equation represented in the (time- 
dependent) Slater determinant basis. The nonlinear 
equation ( 12 ) will be called the orbital equation. In it 



we introduced the one- and two-particle density matrices 
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as well as the mean-field potential operator g r 
coordinate representation reads 



dr . 



1 



(13) 
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which in 
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Particularly in the literature on (single-determinant) 
Hartree-Fock, g rs (r) is also called Coulomb- (r = s) and 
exchange- integrals (r ^ s). Further, in Eq. (12) the elim- 



ination of the Lagrange multipliers led to a projection 
operator 



constructed from M time-dependent orbitals of the form 



<f>n(t)) 



(17) 



As M is typically much smaller than Nb, the correspond- 
ing Slater determinant basis is considerably reduced, 
while the space that can be accessed by the single-particle 
basis is practically identical. In this way, it is possible 
that the orbitals match the physical intuition: an orbital 
could for example be composed dominantly of a bound 
state and, at the same time, also contain an ionized frac- 
tion, as it could be induced by an external field. Orbitals 
adapted in such a way are then likely to yield a more com- 
pact and appropriate Slater determinant basis, so that 
the wavefunction can hopefully be expressed adequately 
by a small set of optimized determinants and not by a 
large static basis. Lastly, the efficiency of the MCTDHF 
approach depends on the degree of correlation present in 
the wavefunction. 



C. Discretization of the orbital equation 



Once the orbitals are discretized through the expan- 
sion in an orthonormal basis (17), upon insertion into 
the orbital equation and projection onto ( ?pk | we readily 



obtain an equation governing the time-dependence of the 
coefficients, 



A I 



P = 1 - E \ <t>™){<t>v 



(16) 



which projects on the orthogonal complement of the sub- 
space spanned by the orbitals. 

Before we proceed in the numerical solution of the 
MCTDHF equations, let us take a comparison of the 
present method with the (full) Configuration Interaction 
(CI) method. The latter also employs a linear expan- 
sion of the wavefunction, however in a basis of static 
Slater determinants. They are constructed from a time- 
independent single-particle basis {| ip^ a )} of size JV& (in- 
stead of M), leading to a size of the determinant basis 
of (jv b )(A^)- As the simulation of ionization processes 
typically requires a large single-particle basis in order to 
provide an adequate description of the scattering states 
(in this work we use up to Nj, ~ 10000), the many-body 
Hilbert space quickly becomes prohibitively large causing 
CI to be applicable at most to two-particle systems. In 
contrast, MCTDHF may provide a much more compact 
representation of the wavefunction by using determinants 
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The superscripts label the number of performed trans- 
formation steps. Let the electron integrals in the time- 
independent basis {| tpk )} be denoted by \) and g, respec- 
tively, and defined analogously to Eqs. ^ and Q. The 
one-body matrix elements are then given by 



hmn (^) 
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FIG. 1. (Color online) Illustration of the finite-element dis- 
crete variable representation basis used to discretize the radial 
part of the wave function. Shown are three equally spaced fi- 
nite elements with 8, 6, 4 Gauss-Lobatto points, respectively, 
(red) dots. Note that the basisfunctions are zero at all grid- 
points except one, and that the functions corresponding to 
the first and last gridpoint have been removed in order to 
satisfy the boundary conditions. 



while the two-body quantities read 
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The quantities without superscripts, Eq. (20) and (24 1, 
are the one- and two-electron integrals in the time- 
dependent basis already encountered in the Hamilto- 
nian The transformations (19)-(25) need to be per- 
formed at each propagation step and constitute the main 
time-consuming part with respect to the single-particle 
basis. In particular, the major effort is the calculation 
of the mean-field potential g rs (r), Eq. (15), which cor- 
responds to evaluation of the first two summations in 



Eq. (23). It grows as 0(N b ) for a general basis, like 



Gaussian- or Slater type orbitals, Sturmians, etc., and 
even though there exist techniques to reduce the effort by 
one order to O(N^) - for instance by using a Cholesky 
decomposition of the two-electron matrix |45| - we found 
that the size of the basis is still limited to, say, Nb ~ 150, 
which is too small for our purpose. 



D. DVR product basis 

In the following, we show how to overcome the afore 
mentioned problems by using an appropriate single- 
particle basis. For the treatment of three-dimensional 



atoms, a spherical expansion provides a description that 
can be much more compact than in cartesian coordinates 
[33] , owing to the fact that the wavefunction is often very 
smooth perpendicular to the radial direction. Hence, we 
use an expansion in spherical harmonics Y/ m , 



(26) 



hit, 



The radial coordinate is described by a discrete variable 
representation (DVR) basis Xfc(r). As the DVR is in 
common use for roughly 30 years in quantum mechanics 
[4"6l |4"7] and was recently successfully used also in many- 
body calculations [48, 49], we refer the reader to the cited 
publications for a general introduction. We only note its 
main advantage, namely the diagonal representation of 
spatially local operators, 



( Xrn \f(r)\ Xn ) = S rnn f(r m ) 



(27) 



where the diagonal entries are simply given by the func- 
tion values on a DVR grid. The DVR thus resembles 
the representation on a finite difference grid, but with 
an important advantage: while finite difference schemes 
converge polynomially with the number of gridpoints, the 
DVR may provide an exponential convergence behaviour. 
This is related to the representation of nonlocal opera- 
tors, like derivative operators, which in general are given 
by full matrices. 

To further take advantage of this feature, but avoid 
full matrices, we particularly use the finite element DVR 
(FEDVR) of Rescigno et al. [SSHS2]- The correspond- 
ing grid of size N ra d consists of finite elements of Gauss- 
Lobatto nodes, see Fig. [I] The advantage in using finite 
elements is that the kinetic energy matrix attains a block 
structure, as matrix elements of one-particle operators 
are nonzero only within adjacent elements. As it is shown 
in [53], by choosing a reasonable distribution and size of 
the finite elements one can avoid large energy eigenvalues 
of the single-particle Hamiltonian, which would make the 
MCTDHF equations more stiff and slow down the time- 
evolution. Further, the positioning of the finite elements 
can be matched to the physical problem, e.g. by employ- 
ing a denser grid near the atom core. 

Having introduced the product basis, the main idea 
to improve the efficiency is now to employ the Poisson 
equation for the calculation of the mean-field potential, 
as it is done routinely in (single-determinant) Hartree- 
Fock implementations, see e.g. [20l|54]. Upon application 
of the Laplace operator to Eq. (15), we readily obtain 



Ag r 



-47r#(r) s (r). 



(28) 



We solve this equation by expansion of g rs ( r ) m the ba- 
sis (26). Note that here the grid-like DVR treatment is 



essential, as it provides sufficient flexibility in represent- 
ing the mean-field potential. The discretized version is 
treated with a method similar to the one used in Ref. |55j . 
It involves the solution of linear systems having a banded 
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structure, which is caused by the finite element basis. 



Thereafter, the integral (22) is calculated using the DVR 



quadrature rule and analytical formulas for the product 
of three spherical harmonics [55]. All in all, we obtain 
a scaling of 0{N ra d N ang ), [here N ang denotes the num- 
ber of angular basis functions] , that is linear with respect 
to N ra d , thus allowing for the application of rather large 
radial grids. We finally note, that particularly for the 
chosen FEDVR + spherical harmonic basis, McCurdy et 
al. have derived a similar method for the calculation of 
the Coulomb repulsion integrals |51j . which may also be 
exploited in the solution of Eq. ( 28 ) . 



E. Numerical implementation 

We briefly mention the basic ideas and properties of 
our program Kiel MCTDHF. It is able to treat both 
fermions and bosons in a unified way where, for the lat- 
ter, an expansion in permanents instead of Slater de- 
terminants is used. For fermions, it provides the spin- 
restricted treatment employed in this work, as well as the 
spin-unrestricted version of the MCTDHF equations. Be- 
side the single-particle basis set used here, several others 
such as Gaussian- and Slater type orbitals, a spheroidal 
DVR grid and some one-dimensional DVRs have been 
implemented, which allow for a treatment of one-, two-, 
and three-dimensional systems. The structure of the 
MCTDHF equations thereby allows for an implementa- 
tion of the different basis sets independent of the wave- 
function part, i.e. of the solution of Eq. (11). The time- 



propagation of the equations of motion is performed with 
the general purpose integrators of Ref. [57 , of which 
we apply most of the times an eight-order Runge-Kutta 
method with adaptive stepsizes. The initial groundstate 
is found via propagation in imaginary time starting from 
a random state. 

For the evaluation of the rhs. of the wavefunction 
equation (11), we mainly follow the Configuration In- 
teraction techniques described in Ref. [32]: the Hamil- 
tonian matrix is calculated using the minimal operation 
count (MOC) method, which utilizes a separation of the 
Slater determinants in products of a- and /3-spin strings. 
Due to the absence of spin-interactions in the Hamilto- 
nian Q, we use only the configurations with vanishing 
spin-projection S z = in the expansion of the wavefunc- 
tion 

As the Hamiltonian and the density matrices have to 
be constructed in each time-step, the matrix elements of 
the excitation operators, Eqs. Q and ([5]), are required 
repeatedly. They are determined at the beginning of 
the run by using a graphical representation method and 
stored in the following. Due to the representation of the 
determinants in spin-strings, only little memory is needed 
for their storage. 

The solution of the orbital equation is found as dis- 
cussed above. Its effort scales as 0(M 2 Ni,N ang ), what is 
much larger than the size of the determinant basis given 



by rf) (* f ) = M 2 - The overall performance of the calcu- 
lations presented here is thus determined by the orbital 
equation. For larger atoms, however, this situation is in- 
verted due to the exponential scaling of the many-body 
basis size. 



F. Time-dependent Configuration Interaction 

As noted above, the helium atom has been subject 
of many investigations based on a direct solution of the 
Schrodinger equation and hence, there are many results 
available which could serve as a reference. Nevertheless, 
we implemented our own direct solution method, in order 
to apply the same single-particle basis as for the MCT- 
DHF calculations and to focus exclusively on the effects 
of the many-body basis. This way, we achieve a more 
consistent comparison, which is particularly independent 
of the single-particle basis and which can be performed 
strictly on a numerical level. 

The actual solution of the Schrodinger equation is 
found using the time-dependent (full) Configuration In- 
teraction (TDCI) method, which has been briefly dis- 
cussed above. Our approach offers the appealing ad- 
vantage, that the basic routines can be adopted from 
MCTDHF without any changes. The time-independent 
eigenvalue equation is solved by a restarted Lanczos algo- 
rithm with full rc-orthogonalization [58] . Subsequently, 
the groundstate is propagated using the short-iterative 
Lanczos method, which basically employs the same al- 
gorithm to approximate the exponential time-evolution 
operator [59] . In a spherical harmonic basis like ( 26 1 , the 



method becomes equivalent to a time-dependent close- 
coupling (TDCC) treatment, the only difference is that 
the angular part is represented by the uncoupled quan- 
tum numbers (h,hi 1^1,^2) instead of the coupled set 
(L,M,h,l 2 ). 



III. NUMERICAL RESULTS 
A. Extraction of physical quantities 

The route towards observables is, in principle, straight- 
forward: given the final wavefunction, all information 
can be obtained by projection on subsets of the complete 
spectrum of atomic eigenstates. For example, single and 
double ionization events can be identified by projecting 
on the corresponding singly and doubly ionized Helium 
eigenstates. The determination of these states, however, 
would require the full diagonalization of the Hamiltonian, 
which is far beyond reach. One must, therefore, resort 
to a projection on approximate states; these include un- 
corrected eigenstates. e.g. antisymmetrized products of 
single-particle eigenf unctions [HI [60] as well as, less fre- 
quently, states that contain a certain amount of corre- 
lation [5J [5]. It is important to note that the results 
from both approaches differ significantly, and that the 
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discussion on the correct treatment has not settled yet 
[3 [5j [UJ . As we are mainly interested in the capabili- 
ties of the MCTDHF scheme, in this work, we apply the 
simpler alternative and neglect correlation in the ionized 
eigenstates. 

In the MCTDHF method, rather than employing the 
many-body wavefunction, it is convenient to use the 
single- and two-particle density matrices which have to be 
supplied at each propagation step in the solution of the 
orbital equation ( 12 ). With them, the expectation values 
of all single- and two-particle operators are available. In 
coordinate representation the corresponding single- and 
two-particle densities are given by [H] 

D(r, t) =J2 D pS) <^( r ' *)^( r < *) . ( 29 ) 

pq 

d(r,r\t) = l - 5^d p9P ,(t)^;(r,t)0,(r,t)^(r',t)^(r' ) t). 

pqrs 

(30) 

Note that for two-particle systems, the two-particle den- 
sity is equal to the absolute square of the wavefunction 
|\&(r, r', t)\ 2 . With the densities, we have access to the 
single- and double ionization yield, defined as the fraction 
of the densities outside a chosen radius Rq, 



Pi(t) 



I 



D(r,t) dr, 



r>R 



P2(t) 



d(r,r',t) drdr' . 



(31) 



(32) 
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Analogous definitions have been used in Ref. |15) . Sim- 
ilarly, it is also possible to study angular distributions 
by integrating only over a specified solid angle. As men- 
tioned above, the formulas ( 31|32 ) can also be viewed as 
resulting from a projection, in which the singly and dou- 
bly ionized continuum states are assumed to be localized 
in different spatial regions. Obviously, this is a crude ap- 
proximation: for instance, bound Rydberg states, which 
extend over a large spatial region, may contribute to the 
ionization yield, and vice versa, delocalized continuum 
states contribute to the bound fraction. Nevertheless, it 
has been possible to relate the thus obtained results to 
experimental data [15] . It is further advantageous, that 
Eqs. ( 31|32 ) are directly applicable also to larger atoms, 
where projection schemes would require a specific set of 
scattering states. 



B. Groundstate 

In this work, the time-evolution is started from the 
groundstate, which itself is found by imaginary time 
propagation (ITP) of a random state according to 
the MCTDHF equations. As the corresponding time- 
evolution operator is non-unitary, at each propagation 
step the wavefunction coefficients have to be renormal- 
ized. In addition, we orthonormalize the orbital expan- 



approximation 


energy [a.u.] 


correlation 


HF 


-2.86178 


0% 


M = 2 


-2.87805 


40% 


M = 3 


-2.88484 


56% 


M = 4 


-2.89135 


72% 


M = 5 


-2.89775 


87% 


M — 7 


-2.89905 


91% 


M = 9 


-2.90017 


93% 


M = 15 


-2.90210 


98% 


CI 


-2.90285 


100% 


exact [ST] 


-2.90372 





TABLE I. Helium groundstate energies for different numbers 
M of time-dependent orbitals, and percentage of the covered 
correlation energy. The difference between our Configuration 
Interaction result and the exact nonrelativistic groundstate 
energy is due to the angular basis, which includes partial 
waves only up to / = 2. 



sion coefficients, although this is in principle not neces- 
sary, owing to the constraint in the action functional ( 10 1 



However, it yields a more stable scheme and prevents 
small numerical errors from growing during the ITP. 

In Tab. [TJ we present the groundstate energies for dif- 
ferent numbers M of MCTDHF orbitals. In each calcu- 
lation, we use a single-particle basis including all partial 
waves up to I — 2, and a radial part which is described 
by three finite elements of length 4 a.u. with 11 Gauss- 
Lobatto nodes, that span the region up to r max — 12 a.u. 
This is a typical gridpoint distribution also employed in 
the time-dependent calculations presented below, i.e. not 
optimized for groundstate calculations. The results are 
compared to the full CI energy, which marks the fully 
correlated result for the chosen single-particle basis. In 
the second row of the table, we give the fraction of the 
correlation energy, defined as the difference between the 
CI and the HF result, that is covered by the respective 
approximation. 

With increasing numbers of self-consistently deter- 
mined MCTDHF orbitals, the energies obviously con- 
verge against the CI reference result. The first correction 
to Hartree-Fock, M — 2, is thereby able to account for 
40 % of the total correlation energy. On the other hand, 
the most accurate approximation considered in this work, 
M = 15, covers 98 % of the correlation energy and dif- 
fers from the CI energy only by 0.8 mHa. This meets the 
(loosely specified) requirements of chemical accuracy, and 
should be appropriate for most applications. 

The observed convergence is substantially slower than 
for one-dimensional model atoms (39J [40] ; we have ver- 
ified by using different initial states and basis sets that 
the groundstate has indeed been reached. Next, we study 
whether the description of time-dependent processes is 
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FIG. 2. (Color online) (a) Single- and (b) double ionization 
yield vs. intensity for a photon energy of 45 eV. Shown are 
the results for different MCTDHF approximations and TDCI 
as well as the results of Nikolopoulos and Lambropoulos, Ref. 



similarly accurate. 



FIG. 3. (Color online) (a) Single- and (b) double total ioniza- 
tion cross sections vs. photon energy, at a fixed field intensity 
of 10 13 W/cm 2 . Shown are the results of selected MCTDHF 
approximations and from TDCI in the same basis, the results 
of Feist et al. as well as the results of experiments by 

Hasegawa et al. |12j and Sorokin et al. [13] (the latter with 
errorbars) . 



C. Two-photon ionization 

In the following we present the results for the two- 
photon ionization induced by electromagnetic fields with 
frequencies in the range 40 to 54 eV, and intensities from 
10 13 to 10 16 W/cm 2 . The fields are linearly polarized 
and have a squared-sine envelope, 



E(t) = E sin(u;£)sin 2 (7rf/T) 



(33) 



for < t < T and zero else. Throughout this work, we 
use a duration T = 100 a.u. (2.4 fs), which corresponds 
to roughly 24 to 32 cycles for the applied frequencies. 
The single-particle basis is discretized on a radial grid 
that extents up to r max — 100 a.u., with finite elements 
of length 4 a.u. containing 11 FEDVR functions each. 
The maximal angular quantum numbers are I = 3 and 



m = 1. For this choice of parameters, we obtain well 
converged results; a change to r max = 150 a.u., I — 4 and 
rn = 2 yielded virtually the same results. Observables 
are extracted at the end of the pulses. The ionization 
radius in Eqs. ( 31|32 ) is chosen to be R = 20 a.u. We 
show results obtained in the length gauge; the velocity 
gauge was found to yield virtually identical results for 
the various test cases we considered. 

In Fig. [2] (a), we plot the single- ionization yields ver- 
sus intensity for a photon energy of 45 eV. At intensities 
lower than 5 • 10 14 W/cm 2 , where perturbation theory 
is applicable, all depicted curves are given by straight 
lines with a slope of one. Over the whole range of in- 
tensities, the TDCI results are well reproduced already 
with the first correction to TDHF, M = 2, while M = 5 
yields full agreement. Fig. [2] (b) depicts the correspond- 
ing double ionization results. They are compared to the 
results of Nikolopolous et al. [5] , which beside a system- 




photon energy [eV] ionization radius R [a.u. ] 



FIG. 4. (Color online) Convergence behaviour of the MCT- 
DHF double ionization cross sections. The successive approxi- 
mations M — 2, 3, M — 6,7 and M — 8, 9 give nearly identical 
results. 



atic shift agree well with our TDCI curve. The TDHF 
result shows a similar qualitative behaviour, but is more 
than one order of magnitude larger. The plotted MCT- 
DHF results deviate from the expected straight lines. 
The M — 2 approximation lies roughly in between TDHF 
and TDCI, and approaches the latter for intensities larger 
than 5-10 14 W/cm 2 . The M = 10 result is for most of the 
considered intensities more accurate; however, it shows 
an unsatisfactory behaviour, particularly at the intensity 
2 • 10 14 W/cm 2 , where a further increase of the density 
leads to an unphysical decrease of the ionization yield. 

The largely overestimated double ionization yield ob- 
tained from TDHF is caused by the nature of the sin- 
gle determinant approximation. Upon insertion of the 
Hartree-Fock two-particle density, which is for two parti- 
cles given by o? HF (r, r') = £>(r)_D(r')/4, in the definition 
(32 1, one obtains Pn = /4. This shows that the double 



ionization is directly related to the single ionization, and 
that, as a consequence, the TDHF approximation fails to 
model the physical picture of double ionization correctly. 
In a similar way, low-order MCTDHF approximations 
should be affected by this mechanism, although with in- 
creasing number of orbitals M, the one- and two-body 
density matrices and, thus, also the ionization yields be- 
come more independent from one another. 

In Fig. |3j we show the total cross sections for the two- 
photon single and double ionization obtained at a field in- 
tensity of 10 13 W/cm 2 . For the used squared-sine pulses, 
the cross sections can be related to the ionization yields 
( |3l]32| through [5] 



CTi [cm 2 ] = 1.032 • 10 
CT 2 [cm 4 s] = 2.28 • 10" 



_4 " 2 Pi 
nl 

uj 3 P 2 
nli 



(34) 



(35) 



FIG. 5. (Color online) Single and double ionization yields for 
selected MCTDHF approximations plotted against the radius 
Ro, at which a particle is considered ionized. The line at 
Ro = 20 a.u. marks the value used in this work. The pulse 
parameters are to — 45 eV and I = 10 13 W/cm 2 . 



is the frequency in eV and Iq is the intensity in W/cm 2 . 
The single-ionization cross section in Fig. [3] (a) shows a 
monotonic decrease with growing photon energy which 
is qualitatively well reproduced by all depicted curves. 
For M = 4, MCDTHF shows a good agreement with the 
TDCI result. In Fig. |] (b) we plot the double ioniza- 
tion cross section. First et us focus on the TDCI results, 
which is compared to the results of Feist et al. [TT] , that 
likewise were obtained via a direct solution of the TDSE. 
The deviation between the curves is caused by the differ- 
ent methods to extract the double ionization yield (see 
for instance Fig. 8 in that reference, where several other 
methods show similar discrepancies). In |llj it has also 



been observed, that the expression ( 32 1 overestimates the 



where n represents the number of cycles in the pulse, w 



double ionization as compared to the one obtained via 
projection on Coulomb waves by roughly 25% (they used 
an ionization radius of R — 70 a.u.); as can be seen from 
our calculations, this ratio is even larger in the case of 
the smaller radius Rq = 20 a.u. used here, particularly 
at frequencies u> > 48 eV. We also depict the results of 
two experiments [121 113] , which are reasonably close to 
our TDCI result. The M = 15 approximation fails to 
give a quantitatively valid description: at small photon 
energies, it underestimates the TDCI results, while for 
photon energies larger than 48 eV the yields become too 
small. Further, the monotonic growth of the TDCI result 
is not reproduced. However, the approximation is able to 
qualitatively describe the increase of the curve near the 
threshold. 

In Fig. |4j we take a closer look on the convergence 
behaviour of the MCTDHF approximations. One can 
see, that TDHF again yields a far too large result for 
small photon energies. The approximations M = 2, 3 and 
M = 6, 7 improve the order of magnitude, but are not 
able to reproduce the increase near the threshold. The 
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FIG. 6. (Color online) Snapshots of the radial density for 
different times t from 100 to 400 a.u. after application of a 
42 eV pulse of intensity 10 13 W/cm 2 , calculated on a grid 
that extends up to r max = 500 a.u. One clearly notices 
a wavepacket leaving the atom, which moves slightly faster 
within the M = 3 approximation than within TDHF. 



approximation M = 8 is the first to yield a reversal of 
this trend around 52 eV. It is further interesting, that the 
pairs of successive approximations M — 2,3, M = 6, 7 
and M = 8, 9 lead to almost identical results. 

A shortcoming of our definition of ionization, 
Eqs. ( 31|32[ ), is the arbitrariness of the ionization radius 
Rq. In Fig. [5] we thus examine its influence on the sin- 
gle and double ionization yields induced by a pulse with 
u) = 45 eV and / = 10 13 W/cm 2 . One can see that for 
all depicted MCTDHF approximations, there is a region 
between r — 8 and r — 40 a.u., in which the yields are 
almost constant. From this, one can expect that for the 
chosen grid dimension and excitation scenario our choice 
of Rq = 20 a.u. is a reasonable one. 

The preceding results were obtained on a grid extend- 
ing to r max = 100 a.u. in order to keep the effort of 
TDCI manageable. As stated before, through the use of 
the Poisson equation, MCTDHF provides a linear scaling 
with respect to the number of gridpoints. For our final 
result, we exploit this fact by using a five times larger 
grid with r max = 500 a.u. with the same finite element 
distribution as before. In Fig. [6j we show snapshots of 
the radial density induced by a pulse with w = 42 eV 
and I = 10 13 W/cm 2 , for different times up to t = 400 
a.u. (~ 10 fs). One clearly notices a wavepacket mov- 
ing away from the atom; at t = 400 a.u. a significant 
part of it starts to get reflected at the grid boundary, 
and subsequently interferes with its outgoing part. Such 
reflections could be avoided relatively easily by applying 
exterior complex scaling or a complex absorbing poten- 
tial. However, this is not necessary for the results pre- 
sented in this paper since, at the time t — 100 a.u., most 
of the wavepacket is still in the region inside 100 a.u. 

One can further observe, that the velocity of the 



wavepacket, which is determined by difference of the pho- 
ton energy and the ionization potential, is larger in the 
M = 3 approximation than in TDHF. This can be ex- 
plained with the well-known fact, that the Hartree-Fock 
approximation - or more precisely the Koopman's theo- 
rem - overestimates the exact value of the ionization po- 
tential. The MCTDHF approximation improves on this, 
and the resulting smaller ionization potential leads to a 
larger electron excess energy. 



IV. DISCUSSION 

The numerical investigation of the two-photon ioniza- 
tion of Helium presented in this paper revealed much 
of the capabilities and limitations of the MCTDHF for- 
malism. Through comparison with time-dependent (full) 
Configuration Interaction calculations using the same 
single-particle basis, we were able to critically assess the 
accuracy of different MCTDHF approximations on the 
many-body level. The TDCI results itself were com- 
pared to results from the literature and showed a reason- 
able agreement. The conclusions drawn in the following 
are likely to apply also to other excitation scenarios and, 
more importantly, to atoms with more electrons, which 
is the desired area of application of MCTDHF. 

First, we observed that the single-particle ionization 
is very well reproduced already with a small number of 
time-dependent orbitals (M ~ 4), whereas TDHF yields 
qualitative agreement. The effect of MCTDHF is thereby 
to provide a correlated background on which the ioniza- 
tion process takes place and which - in contrast to single- 
active electron calculations - is self-consistently found 
and adapted during the temporal evolution. 

The double ionization is - for the considered photon 
energies - a highly correlated process, and thus a severe 
test for MCTDHF. Correspondingly, its correct descrip- 
tion requires a large determinant basis. For the present 
case of Helium, the maximum number of 225 Slater deter- 
minants (M = 15) is not sufficient to achieve converged 
results, yet it suffices for a qualitative agreement with 
the fully correlated TDCI result (for which more than 6 
million Slater determinants were used). 

One surprising point we observed is that the MCT- 
DHF results become better with increasing intensity, a 
behaviour totally different from approaches based on per- 
turbation theory, such as LOPT. Accordingly, the double 
ionization is described better if the ionized fraction con- 
stitutes a significant part of the total particle number. 
On the contrary, the time-dependent variational princi- 
ple performs worse in modelling a small ionized fraction, 
in favour of giving a well description in the vicinity of 
the atom. For atoms with more electrons, and for a given 
number M of time-dependent orbitals, we however expect 
an improved accuracy of the double ionization as com- 
pared to the Helium case, due to the largely increased 
Slater determinant basis, while the difficulties encoun- 
tered with low-order corrections are hopefully shifted to- 
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wards higher-particle ionization yields. 

The present numerical study of the MCTDHF method 
focused on the two-photon ionization of Helium since here 
the time-dependent Schrodinger equation is directly solv- 
able and direct comparisons are possible. Based on these 
results, our future work will concentrate on the applica- 
tion of the MCTDHF method to the photoionization of 
larger atoms and molecules where direct solutions of the 
Schrodinger equation are not (yet) feasible. 



respect to the orbitals to vanish: 
S 



= 



{Cn},{|<M} 




(Al) 



^2 d pi rs ( <M 5 1 &i ) | <t> q ) r ~ ^?(*) 1 0« ) • 

qrs ^' ^ ^ ) q 
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Appendix A: Spin-restricted MCTDHF equations 

The following derivation of the MCTDHF equations 
is mainly inspired by the work of Alon et al. [33], who 
derived equations for fermions in a spin-unrestricted fash- 
ion where the spin was treated implicitly, i.e. absorbed 
in the orbital index. Here, we present a spin-restricted 
formulation. The appearance of the equations and the 
mathematical procedure is, apart from notational issues, 
essentially identical. Yet, there are differences in the in- 
terpretation: in the following spin-restricted formulation, 
2M spin-orbitals are described by M spatial orbitals, 
\<t> P a{t)) = \<j) p p{t)) =: \4> p (t)). The indices of the 
matrix elements and summations correspondingly range 
from 1 to M. The unrestricted treatment in Ref. 33 is 
formulated in terms of the M a -\-Mg spin-orbitals. There, 
the indices run from 1 to M a +Mp and roughly one half of 
the matrix elements are equal to zero, due to the orthog- 
onality of a and ft spin-functions. This is however is not 
indicated in the reference, as the spin summation in the 
electron integrals is lacking. A further difference between 
the two approaches is in the appearance of the density 
matrices, Eqs. (13) and (14 1, which arc here composed 



of summations over the internal spin variables, whereas 
in Ref. |33j they are set up in terms of the elementary 
spin-orbital excitation operators. 

Let us turn to the actual derivation, and first to the or- 



Applying the projector | cj> m )( <f> m |, summing over m and 
solving for the term with the Lagrange multiplier /x pm , 
we obtain: 



^AVn(i)|0rn) = ^ | <f>m ) ( <§>m \ x 
d 



(A2) 



h 



'dt 



) ^ dpq rs Q rs | (fig ^ / 
qrs ) 



In the two Eqs. (Al) and (A2), the curly brackets are 



equal. Insertion of the latter in the former thus yields 
d 



h — i- 



dt 



f ->q ) + £ d pqrs g rs | 4> q ) > , 

qrs ) 



(A3) 



where P was defined in Eq. (16). Solving for the term 
with the time derivative, multiplying by (D" 1 ),^ and 
summing over p yields: 



?P 



4>n ) — P < h | 4>n ) + £ (D l ) np d pqrs g rs | (j)q ) > 
I pqrs ) 



(A4) 



In order to obtain explicit equations, we apply a unitary 
transformation among the orbitals after which 



d 

dt 



= o 



(A5) 



is satisfied. This causes the projector on the lhs. to van- 



ish and leads immediately to the orbital equation ( 12 ) 



For the coefficient equation we set the derivative with 
respect to the coefficients equal to zero. Therefore, we 
insert the MCTDHF expansion |9]) into the action func- 
tional to obtain a form which explicitly depends on the 
coefficients. The derivative then directly yields 



= s[{C m },{|<M} 



S 

.dC n 



dt 



i— m 

dt 1 ' 



(A6) 



bital equation. In the action functional (|10|), we express With the property jA5l), also the Slater determinant ma- 



the expectation value of H — i^, in terms of the den- 



sity matrices and require the functional derivative with 



trix element of the time derivative vanishes, and we ob- 
tain the wavefunction equation. 
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Practically, the constraint ( A5 1 avoids rotations within equations 



the subspace spanned by the orbitals. If such a rotation 
is to be performed, only the wavefunction coefficients are 
affected, while the single-particle subspace is varied only 
if unavoidable. However, the applied unitary transforma- 
tion is only a special representative of a broader class of 
constraints which read 



d_ 
dt 



Umn (^) ; 



(A7) 



where U mn (t) is an arbitrary time-dependent hermitian 
matrix. This leads to a slightly different form of the 
MCTDHF equations, see Refs. [21133]. For instance, 



by choosing U mn (t) 



(t) one obtains the working 



iC n (t) = ^2 ( n | - ^2 9 P qrs e pqrs | m) C m (t) , (A8) 



i\4> n ) = h(t)\4> n ) + P 



{ E (p- 1 ) 

^ pqrs 



t dpqrs 9rs \ } 



(A9) 



We tested them as an alternative for the real-time prop- 
agation and observed a roughly ten percent slower per- 
formance, in contrast to what is reported for MCTDH 
calculations in Ref. 28J. This is likely to be caused by 
the dominant contribution of the solution of the orbital 
equation to the total effort. 

We finally note that the previous derivation relies heav- 
ily on the orthonormality of the underlying MCTDHF 
orbitals, which is enforced by the Lagrange multipliers in 
the action functional ( 10 ). If we allowed for nonorthonor- 



mal orbitals, several complications would arise, as can be 
seen, for instance, in the derivation of the unrestricted 
bosonic Hartree-Fock approximation [62 . Particularly, 
and most intriguingly, the derivative of the density ma- 
trices with respect to the orbitals would no longer vanish. 
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